This tutorial is part of FEMlium.

import numpy as np
import pyproj
import dolfinx
import dolfinx.plot
import pyvista
import folium
import femlium
from utils import gmsh_to_fenicsx

Auxiliary function to get a folium Map close to Lake Garda.

def get_garda_geo_map(boundary_icons=False):
    # Add map close to Lake Garda
    geo_map = folium.Map(location=[45.6389113, 10.7521368], zoom_start=10.3)

    # Add markers
    if boundary_icons:
        location_markers = {
            "Sarca": [45.87395405, 10.87087005],
            "Mincio": [45.43259035, 10.7007715]
        }
        location_colors = {
            "Sarca": "red",
            "Mincio": "green"
        }

        for key in location_markers.keys():
            folium.Marker(
                location=location_markers[key],
                tooltip=key,
                icon=folium.Icon(color=location_colors[key])
            ).add_to(geo_map)

    # Return folium map
    return geo_map

get_garda_geo_map()

Make this Notebook Trusted to load map: File -> Trust Notebook

Read the mesh, the subdomain markers and the boundary markers from file with dolfinx.

mesh, subdomains, boundaries = gmsh_to_fenicsx("data/garda.msh")

Plot the mesh using pyvista.

def dolfinx_to_pyvista_mesh(mesh):
    num_cells = mesh.topology.index_map(mesh.topology.dim).size_local
    cell_entities = np.arange(num_cells, dtype=np.int32)
    pyvista_cells, cell_types = dolfinx.plot.create_vtk_topology(mesh, mesh.topology.dim, cell_entities)
    grid = pyvista.UnstructuredGrid(pyvista_cells, cell_types, mesh.geometry.x)
    return grid

def pyvista_mesh_plot(mesh):
    grid = dolfinx_to_pyvista_mesh(mesh)
    plotter = pyvista.PlotterITK()
    plotter.add_mesh(grid)
    plotter.show()

pyvista_mesh_plot(mesh)

Plot the subdomain markers using pyvista.

def pyvista_subdomains_plot(mesh, subdomains):
    grid = dolfinx_to_pyvista_mesh(mesh)
    grid.cell_arrays["Marker"] = subdomains.values
    grid.set_active_scalars("Marker")
    plotter = pyvista.PlotterITK()
    plotter.add_mesh(grid)
    plotter.show()

pyvista_subdomains_plot(mesh, subdomains)

WARNING:py.warnings:Use of `cell_arrays` is deprecated. Use `cell_data` instead.

Define a pyproj Transformer to map between different reference systems, because the points read from file are stored a $(x, y)$ pairs in the EPSG32632 reference system, while the map produced by folium is based on (latitude, longitude) pairs in the EPSG4326 reference system.

transformer = pyproj.Transformer.from_crs("epsg:32632", "epsg:4326", always_xy=True)

We define a mesh plotter for meshes in dolfinx format, which is implemented in femlium.DolfinxPlotter.

dolfinx_plotter = femlium.DolfinxPlotter(transformer)

We use the dolfinx_plotter to draw the mesh on top of the geographic map.

geo_map = get_garda_geo_map()
dolfinx_plotter.add_mesh_to(geo_map, mesh)
geo_map

Make this Notebook Trusted to load map: File -> Trust Notebook

We may change the color and the weight of the line.

geo_map = get_garda_geo_map()
dolfinx_plotter.add_mesh_to(geo_map, mesh, face_colors="red", face_weights=2)
geo_map

Make this Notebook Trusted to load map: File -> Trust Notebook

Furthermore, we may set the colors and the weights of the face representation to depend on the markers associated to each segment.

geo_map = get_garda_geo_map(boundary_icons=True)
face_colors = {
    0: "gray",
    1: "blue",
    2: "red",
    3: "green"
}
face_weights = {
    0: 1,
    1: 2,
    2: 5,
    3: 5
}
dolfinx_plotter.add_mesh_to(
    geo_map, mesh, face_mesh_tags=boundaries, face_colors=face_colors, face_weights=face_weights)
geo_map

Make this Notebook Trusted to load map: File -> Trust Notebook

Cells can be colored as well, with a uniform color or depending on the cell markers. We start from a uniform color.

geo_map = get_garda_geo_map()
dolfinx_plotter.add_mesh_to(geo_map, mesh, cell_colors="orange")
geo_map

Make this Notebook Trusted to load map: File -> Trust Notebook

We also show the case of colors being set from cell markers. There are two cell markers in this mesh, equal to 1 for the region close to the shoreline (colored in purple) and 2 for the rest of the domain (colored in yellow).

geo_map = get_garda_geo_map()
cell_colors = {
    1: "purple",
    2: "yellow"
}
dolfinx_plotter.add_mesh_to(geo_map, mesh, cell_mesh_tags=subdomains, cell_colors=cell_colors)
geo_map

Make this Notebook Trusted to load map: File -> Trust Notebook

Once can use colors associated to both cell and face markers on the same plot.

geo_map = get_garda_geo_map(boundary_icons=True)
dolfinx_plotter.add_mesh_to(
    geo_map, mesh,
    cell_mesh_tags=subdomains, face_mesh_tags=boundaries,
    cell_colors=cell_colors, face_colors=face_colors, face_weights=face_weights)
geo_map

Make this Notebook Trusted to load map: File -> Trust Notebook

In order to define a simple scalar field, we compute the centroid of the domain.

centroid = np.mean(mesh.geometry.x[:, :2], axis=0)

We may plot the centroid on top of the mesh.

geo_map = get_garda_geo_map()
dolfinx_plotter.add_mesh_to(geo_map, mesh)
folium.Marker(location=transformer.transform(*centroid)[::-1], tooltip="Centroid").add_to(geo_map)
geo_map

Make this Notebook Trusted to load map: File -> Trust Notebook

The scalar field is defined as $s(\rho, \theta) = \frac{\rho}{\sqrt{1 - 0.5 \cos^2 \theta}}$, and is interpolated on a $\mathbb{P}^2$ finite element space. Here $(\rho, \theta)$ are the polar coordinates centered at the centroid.

scalar_function_space = dolfinx.FunctionSpace(mesh, ("CG", 2))

def scalar_field_eval(x):
    rho = np.sqrt((x[0] - centroid[0])**2 + (x[1] - centroid[1])**2)
    theta = np.arctan2(x[1] - centroid[1], x[0] - centroid[0])
    return rho / np.sqrt(1 - 0.5 * np.cos(theta)**2)

scalar_field = dolfinx.Function(scalar_function_space)
scalar_field.interpolate(scalar_field_eval)

We next show a filled contour plot using pyvista.

def pyvista_scalar_field_plot(mesh, scalar_field):
    grid = dolfinx_to_pyvista_mesh(mesh)
    grid.point_arrays["Scalar field"] = scalar_field.compute_point_values()
    grid.set_active_scalars("Scalar field")
    plotter = pyvista.PlotterITK()
    plotter.add_mesh(grid)
    plotter.show()

pyvista_scalar_field_plot(mesh, scalar_field)

WARNING:py.warnings:Use of `point_arrays` is deprecated. Use `point_data` instead.

In order to plot a field on a geographic map, we use again the dolfinx_plotter. We may plot a filled contour plot on the geographic map.

geo_map = get_garda_geo_map()
dolfinx_plotter.add_scalar_field_to(geo_map, scalar_field, mode="contourf", levels=15, cmap="jet")
geo_map

Make this Notebook Trusted to load map: File -> Trust Notebook

Similarly, we can also use (unfilled) contour plots.

geo_map = get_garda_geo_map()
dolfinx_plotter.add_scalar_field_to(geo_map, scalar_field, mode="contour", levels=15, cmap="jet")
geo_map

Make this Notebook Trusted to load map: File -> Trust Notebook

One may also combine mesh plots and solution plots.

geo_map = get_garda_geo_map()
dolfinx_plotter.add_mesh_to(geo_map, mesh, face_colors="grey")
dolfinx_plotter.add_scalar_field_to(geo_map, scalar_field, mode="contour", levels=15, cmap="jet")
geo_map

Make this Notebook Trusted to load map: File -> Trust Notebook

We next define a vector field $\mathbf{v}(\rho, \theta) = \begin{bmatrix}-\rho \sin \theta\\\rho \cos\theta \end{bmatrix}$.

vector_function_space = dolfinx.VectorFunctionSpace(mesh, ("CG", 2))

def vector_field_eval(x):
    rho = np.sqrt((x[0] - centroid[0])**2 + (x[1] - centroid[1])**2)
    theta = np.arctan2(x[1] - centroid[1], x[0] - centroid[0])
    values = np.zeros((2, x.shape[1]))
    values[0] = - rho * np.sin(theta)
    values[1] = rho * np.cos(theta)
    return values

vector_field = dolfinx.Function(vector_function_space)
vector_field.interpolate(vector_field_eval)

We first see a plot obtained with pyvista, which shows both the magnitude of the vector field and its representation using glyphs.

def pyvista_vector_field_plot(mesh, vector_field):
    grid = dolfinx_to_pyvista_mesh(mesh)
    values = np.zeros((mesh.geometry.x.shape[0], 3))
    values[:, :2] = vector_field.compute_point_values()
    grid.point_arrays["Vector field"] = values
    grid.set_active_vectors("Vector field")
    plotter = pyvista.PlotterITK()
    plotter.add_mesh(grid)
    glyphs = grid.glyph(orient="Vector field", factor=1e-1)
    plotter.add_mesh(glyphs)
    plotter.show()

pyvista_vector_field_plot(mesh, vector_field)

We may obtain contourf or contour plots of the magnitude of the vector field.

geo_map = get_garda_geo_map()
dolfinx_plotter.add_vector_field_to(geo_map, vector_field, mode="contourf", levels=15, cmap="jet")
geo_map

Make this Notebook Trusted to load map: File -> Trust Notebook

geo_map = get_garda_geo_map()
dolfinx_plotter.add_vector_field_to(geo_map, vector_field, mode="contour", levels=15, cmap="jet")
geo_map

Make this Notebook Trusted to load map: File -> Trust Notebook

Also a quiver plot can rendered on top of the geographic map.

geo_map = get_garda_geo_map()
dolfinx_plotter.add_vector_field_to(geo_map, vector_field, mode="quiver", scale=1e-1, cmap="jet")
geo_map

Make this Notebook Trusted to load map: File -> Trust Notebook